General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 




- 1 *jC 02 

Jl.CiJL 

1 JJ 


tbAjA-ci - 1 7 J55t>) in 
Urtt u LLiL, 1 L * t FLi k o) 
>1 L 1 ii U J o^» U i: * v . ) 


U J/ J 4 


*»i- h iSil z j^tntUT 
*■ J V AJ^/rii A U 1 
v-i>U 



InstiU :e for 
Numerical Methods 
in Engineering 





!\]f\Grl^ ■ 

HIGH SPEED INVISCID COMPRESS INFLOW 
BY THE FINITE ELEMENT METHOD 

0. t C. Zienkiewicz, R. Ltihner S K. Morgan 

C/R/ 477/84 May 1 984 


Paper presented at the 5th MAFELAP Conference, 
Brunei University, 1-4 May, 1984. 



HIGH SPEED INVISCID COMPRESSIBLE FLOW 
BY THE FINITE ELEMENT METHOD 

O.C.Zienkiewicz, R.Lohner and K. Morgan 

University of Wales, Swansea, Great Britain 


1. INTRODUCTION 

The numerical solution of compressible flow problems has 
received much attention over the past thiry years due, to a 
large extent, to the interests of the aerospace industry. The 
solutions of such problems are cha racteri sed by the appearance 
of discontinuities, such as shock waves, in the flow field and a 
major topic of attention has been the deveopment of numerical 
techniques which are able to adequately resolve such phenomena. 
A recent paper by Woodward Bnd Collella [1] gives an excellent 
survey of the existing 'state of the art' and compares the 
performance of various widely used algorithms for certain test 
problems. The algorithms considered are finite 
difference/finite volume based and utilise either artificial 
viscosity [2], linear hybridization [3] or expLicit 
nonlinearity [43. All the schemes considered are shown to 
possess certain advantages and disadvantages. 

The finite eLement method has only recently made its 
appearance in this area, but it is expected that it will make a 
significant contribution because of the great geometrical 
flexibility which is inherent in the method. A recent paper by 
Hughes [5] investigates the approach of expLicit nonlinearity 
in a finite element context. The present authors have made some 
initiaL studies [6-3] of high speed inviscid flow problems in 
which they have used the finite eLement method and an explicit 
time stepping algorithm which is based on the Tay Lor-Ga lerkin 
schemes of Donea et aL [9-12], with an appropriate artificial 
viscosity [133. In this paper, we combine this solution 
algorithm with an automatic mesh refinement process which is 
designed to produce accurate steady state solutions to problems 
of inviscid compressible flow in two dimensions. The results 
of two test problems are incLuded which demonstrate the 
excellent performance characteristics of the proposed 
procedures . 


S.-A SINGLE STEP ALGORITHM, 


Inviscid compressible flow of an ideal gas in two dimensions 
is governed by the Euler equations, which can be written in 
vector form as 
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Here p, f> and e denote the pressure, density and specific 
internal Ly energy of the fluid respectively while UL 1 is the 


velocity component 


in 


di rection 


of a Cartesian coordinate 


system. In addition, the summation convention is employed with 
denoting the Kronecker delta. 


A single step algorithm for the solutiion of (2.1) has been 
fully described elsewhere (6 3 and is based upon the 
Tay lor-Ga lerkin methods of Donea et al [9-123. A brief 
descripion of the solution procedure will be given here for the 
sake of completeness. 


Using a Tay Lor series expansion about time t 


u°*'= y^-v-At v 




t gives 

(2.43 


correct to second order, where t” 4- = t* + At and a superscript 
n denotes an evaluation at time t = t n . Eliminating the time 
derivatives via (2.1) leads to the time-stepping scheme 
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The region^- over which the solution is required is discretised 
using linear 3-noded triangular eLements and a Galerkin finite 
element solution procedure is appLied to (2.5) using 
approximations 

u ^ o* = S 

5 * 5 *“ f2 - 7 > 

Here denotes the piecewise linear shape function associated 
with node m and P c the piecewise constant shape function 


associated with eLement e. The resulting matrix equation 
system takes the form 


G* h - O') = f" 


where U T — ( Uj. V (2.9) 

— v - s ~ I / • 

and M is a standard mass matrix. An explicit solution 
procedure for (2.83 is adopted. For problems involving strong 
shocks an artificial viscosity term is included. The form 
used is due to Lapidus Cl 3 3 and replaces the quantities 

calculated from (2.8) by smoothed values 


u£ = ur- 
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where CL. is a constant in the range 0.56^»£2 and is a 


representati ve area for element e. 


3 A TWO STEP ALGORITHM 


Finite difference workers have consistently avoided the use 
of single step explicit algorithms because of the computational 
expense involved in performing matrix multiplications of the 
form required in (2.5). They have favoured instead two-step 
methods which are designed to avoid this requirement. A finite 
eLement two-step algorithm may.be produced which has certain 
features in common with the finite difference scheme of 
Burstein [14]. Alternative finite element two step algorithm 
has recently been proposed by Miner and Skop [153. 
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-3.1 The First Step 

Using a lay Lor series expansion about time t=t n gives 
i , .n .x P 

u ~ y. -v (3.D 

correct to first order, where t* 44 = t^-H^fc, and replacing the 
time derivative from (2.1) produces the expression 


r = o n - 


13.2) 


Employing triangular elements with piecewise Linear and 
piecewise constant shape functions as previously defined, 
approximations are taken in the form. 

u n % u** = 2 Um Hm 


(3.3) 
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These approximations are m substituted into (3.2) and a 
weighted residual statement (16) is formed using the weighting 
function • The result is to give immediately as 
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where Ae. denotes the area of eLement E and denotes the 

surface of this eLement. It may be observed that U ^ can be 
obtained quickly as no assembly of element conributions is 
required to form the right hand side of (3.4). 

3.2 The Second Step 

The second step begins by making a new Taylor expansion and 
writing 
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and again replacing the time derivative leads to 
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With the approximations 
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an appropriated wieghted residual form for (3.6) is then 


where & = (Z|,{ v ) is the unit outward normal to the boundary P 
of $. . Inserting (3.7) into (3.B) then gives 


mco~'- 0-). r - 

A 

where H is defined in (2.9) The solution of this equation 
completes the second step. Again the smoothing of (2.10) is 
applied for problems involving strong shocks. It should be 
noted that for Linear problems where 



and A^is constant, (3.8) combined with (3.9) produces exactly 
the single step equation system (2.8) whenK. -Yl. 

This has been confirmed computationally by applying both 
methods to the solution of pure advecticrs problems involving a 
cone-shaped profile in two dimensions. Similar problems have 
been studied extensively in recent years [9, 17, 183. For this 
probLem, the governing equation is the scalar equation. 

■fcU K TV) 

■fct + 'd’Xfc = O (3.11) 

where (A |r A z ) denotes a specified velocity field and U denotes 
the concentration. Fig. la shows the initial profile and Fig. 1b 
illustrates the numerical solution produced after one complete 
rotation in an incompressible rotating velocity field. A 
problem involving the transportation out of the domain of 
interest of an initial cone in an uniform velocity field was 
also investigated. Fig. 1c shows the profile of the 
concentration as the cone is about to Leave the region 
completely. The quality of these results is excellent and they 







were produced without the artificial viscosity term of (2,10) 

i . e. Cp = 0. 

Preliminary tests show the single step and two step 
algorithms to be producing similar results for practical 
probLems involving the non-linear equation set (2.1) - (2.3), 
with the two step method proving to be computationally more 
efficient . 

4. DOMAIN SPLITTING 

It is well-known that solutions of high speed inviscid 
compressible flow problems exhibit narrow regions of rapid 
change (e.g. shocks) which are embedded in Larger regions where 
the solution is smooth. The economical method of producing an 
accurate solution is to use e discretisation which has many 
small elements in the region of rapid change and larger 

elements elsewhers. However, the small elements used might 
then require that a correspondingly small value of the 
timestep, t, ha rs to be employed to ensure stability of the 
algorithms described in the previous section, whereas the 
Larger elements could tolerate the use of a bigger timestep. 
The remedy adopted by the authors is to split the domain into 
regions in which different timestep sizes can be used. This 
method, which bears certain resemblances to the 

explicit-explicit procedures of Belytschko et al [193 and Liu 
[203 is described in detail elsewhere [73 and will not be 
considered further here. It shouLd be noted however that this 
procedure is performed completely automatically by the computer 
code at prescribed times during the analysis. 

5. ADAPTIVE MESH REFINEMENT 

In general, an analyst will have no a priori knowledge of 
the location of those areas of the solution domain in which 
Large gradients will exist. An ideal computational algorithm 
would then require the ability to automatically refine the 
finite element mesh, in zones of high gradients, as the 
computation proceeds. The geometric flexibility of the linear 
triangular element makes it ideaLLy suited to refinement 
processes of this type. 

Adaptive refinement techniques ore currently receiving much 
attention [21,223 and two general classes may be i denti fi ed .viz 
a priori methods and a posteriori methods. Although many 
applications employing a priori grid refinements have been 
reported, at the present time this approach appears to compare 
unfavourably with a posteriori methods with regard to CPU time 
requirements [213. For this reason, we have concentrated on 
the use of a posteriori techniques. Here the soLution is 


obtained on an original mesh and, at a certain time, the mesh 
.is refined according to soma strategy. The solution then 
proceeds, with further refinements being made as requiud, The 
authors have investigated two alternative refinement stategies 
viz mesh movement and mash enrichment. 

5.1 Hesh Movement 

Here the total number of nodes and eLements remains 
constant, but the location of the nodes is changes in order to 
achieve a better overall distribution, Full details of the 
strategy adopted for moving the mesh and handling some of the 
subsequent problems which may arise have been given in detail 
elsewhere [8], This approach alone does not appear to possess 
the degree of flexibility which would be required for general 
aerodynamic problems and so it has been replaced by mesh 
enrichment . 

5.2 Mesh Enrichment 

In this case the original iwafch is held fixed, while 
hierarchical elements [IB] or simply more elements are added. 
Both strategies have their merits and disadvantages. If 
further degrees of freedom are added hierarchicia l ly , the oLd 
shape functions are trained end the new shape functions are 
retained and the new shape functions are introduced for 
relative values of the unknown. Fig. 2 shows the comparison of 
this strategy with the classical enrichment method of adding 
more eLements. The advBtages of the herarchical technique are 
[i ) certain matrices may be re-used, so avoiding some 
recomputation (ii) multigrid solvers may be implemented 
naturally (iii) no geometrical problems are encountered in the 
transition zone between refined and unrefined portions rf the 
mesh. The disadvantages of the hierarchical approach are 
(i) more stiffness coefficients appear as all the refinement 
levels are interconnected lii) the program structure becomes 
very invoLved and, it is expected that, efficient vectori sati on 
will be virtually impossible. 

The classical way of enriching grids is simply to add more 
elements. In this case, all the advantages of the hierarchical 
method are lost whle the disadvantages disappear. However, 
when a detailed comparison was mado of the two possibilities 
[23], it was decided that classical enrichment would be more 
efficient at this stage for general transient problems and so 
this method wae implemented in the computer code. 

5.3 Error Estimate 

Consider a single scalar equation with unknown U and a 


corresponding finite element approximation U*. The aim of any 
adaptive mesh refinement Is to attempt to minimise the highest 
error occurring in the elements representing the domain i.e. 


£ a II U- 0 * 1 ^ (5.1) 
e K 

should be minimised, where 11 ll^ denotes some suitable norm, 
This criterion leads to the requirement that the error should 
be evenly distributed across the domain, which amounts to the 
requirement 


H U- U*l k £= Cxmsbwt 


(5,2) 


For elliptic problems, the appropriate norm appears naturally 
as the energy norm and a fairly complete and general theory of 
error estimators is available [24,25], However, far less 
research has been directed towards hyperbolic problems and a 
detailed theory is still lacking. Here wb attempt to derive e 
simple, yet reliable, error criterion for hyperbolic problems 
and linear elements. 

It has been shown by Strang and Fix [26] that a natural 
norm for first order hyperbolic equations is the defined by 


I Vi.* - 


(5,3) 


If the solution of hyperbolic equations is viewed as an 
approximation problem along characteristics then It can be 
shown [26] that 

* *1.?* iw»r 




where h e is a representative element legnth, and using the 
norm of (5.3) gives 

* (k. 1 1 u nt 

The condition of (5.2) is then replaced by the requirement that 

he 1 ^ 1 1 ~ C**StA*b (5,6) 

or, since the exact soLution U. is unknown, the practical 
requirement becomes 
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he 11 ~ Co ^Wt (5.7) 

If linear eLements ers being employed it seems naturat to 
choose L » 1 and to define 

- uj («‘W i 

If the eLements ere not badly deformed then £5.7) may be 
approximated by the requirement that 

- [\e \ — \Wv\ 0* V- (5.9) 

or ^ /% for nne, two or three dimensions 


5.4 Identification of Regions to be Refined. 

The first step in any mesh refinement algorithm is the 
identification of the regions of the domain in which refinement 
is to ba employed. For (5.9) it is clear that all elements 
which satisfy 


r 

where f> 
respectively. 




(5. ID) 


should be refined further. The constant S in this equation 
acts as a threshold value for the refinement process. Having 
identified in this fashion the patches of eLements that are to 
ba refined, the boundary nodes of these patches are obtained. 
In order to avoid ’saw-tooth 1 type boundary shapes, any element 
which has alL its node on these boundaries is itueLf included 
among the List of eLements to be refined. 


For systems of equations, a ’key-variable’ has to be 
identified in order to employ (5.10). The choice of this 
variable is not obvious but the density has been chosen for the 
computations reported here. 

5.5 Subdivision of the Elements. 


During the subdivision process, three different types of 

element can appear according to the number of nodes of the 

element which lie on the boundary P between the region which 

s 

be subdivided and the region which is to remain 
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unchanged. If the eLement has no nodes on then the element 

is simply divided into four elements as shown in Fig. 3a. If 
the element has one node on P then the element is subdivided 

into two elements as shown in Fig. 3b. Finally, if the element 
has more than one node on then no subdivision is performed 

as shown in Fig, 3c. When this process has been carried out, the 
new nodal coordinates Bnd the new nodal unknowns are calculated 
while the new element connect i vites are formed. The 
subdivision process is completed by identifying new nodes which 
Lie on the boundary of the physical domain and assigning to 
these nodes the appropriate boundary condition marker. 

5.6 Mesh Quality Enhancement 

In many practical situations the straightforward subdivision 
strategy described above wiLL prove satisfactory but, under 
certain circumstances, badly deformed elements can appear. 
Three simple devices hve thus been implemented in an attempt to 
improve the quality of the/ refined mesh and these are described 
below. The resulting meshes are generally free of the 
geometrical difficulties usually encountered when meshes are 
enriched by st raightforward subdivision. In addition, these 
strategies can easi Ly be extended to three dimensions if linear 
tetrahedral elements are employed. 

5.6.1 Smoothing With the refined grid, the sides of the 

elements are replaced by springs of unit stiffness. For badly 

deformed eLements the resulting nodal forces will not be in 

equilibrium, whereas for regions of well-formed elements the 
resulting nodal forces wilL nearLy vanish. A relaxation 
procedure [8] is adopted which moves the nodes until 

equilibrium of the nodal forces is achieved. Five steps of 
this process are normally employed and this ensures a local 
smoothing of the mesh. 

5.6.2 Avoidance of Successive Subdivision into Two Badly 

deformed elements can appear when an eLement that has been 

initially divided into two elements as in Fig. 3b is again 
subdivided in the same fashion. The easiest method of avoiding 
this second subdivision is to enlarge the region which is to be 
refined by adding to it all elements surrounding the eLement 
under consideration. The eLement under consideration can then 
be subdivided into eight elements, thus enhancing the 

regularity of the mesh. This process is iLLustrated in Fig. 4. 

5.6.3 Removal of Badly Deformed Elements If badly ‘-formed 
elements are still present, after avoiding 'cessive 




subdivision into two and after smoothing the mesh, then these 
elements are removed from the computation. This algorithm is 
detailed elsewhere [83 and need not be repeated here. 

6. RESULTS 

Results will be presented which show the application of the 
methods described in this paper to two separate problems viz. 
steady supersonic flow past a wedge and steady supersonic flow 
past a nose cone. As both problems are steady, the transient 
algorithms described previously are used here as a device for 
integrating to steady state. 

6.1 Supersonic Flow Past a Wedge 

The solution domain and the initial discretisation is shown 
in Fig. 5a. Along AB the flow variables are held fixed at the 
values f = 1.185, u ( = 1080, u L = 0, e = 787350 giving b free 
stream Mach number of three. The boundary ED is a solid wall 
while 'natural 1 conditions are applied along BC, CD and AE. 
The exact solution for this problem consists of an attached 
plane shock at an angle of 37.5° to the 3^-axis with uniform 
flow ahead of and behind the shock. 

The density distribution after 100 global timesteps is shown 
in Fig. 5b, and, aLthough the required features are present, the 
quality of the soLution is not good because of the coarseness 
of the initial discretisation. The mesh is therefore 
automatically refined at this stage, as shown in Fig. 6a, and 
the soLution is advanced in time. After a further 100 global 
steps the density distribution is as illustrated in Fig. 6b and 
the improvement in the quaLity of the soLution is apparent. A 
further refinement was performed, producing the mesh of Fig. 7a, 
whch produced the density distribution shown in Fig. 7b after a 
further 100 global steps. The quality of this solution is 
excellent when compared with the behaviour of the exact 
so luti on . 

6.2 Supersonic Flow Past a Nose Cone 

The solution domain and initial discretisation for this 
example is shown in Fig.B. Along AB the flow variables are 
held fixed at the values f = 1, u t = 1.028, u^ - 0, e = 1 
giving a free stream Mach number of two. The boundary CD is a 
solid wall and 'natural' conditions are applied over BC and DA. 
The initial finite eLement discretisation is shown in Fig. 8 and 
the solution was advanced 100 global timesteps on this mesh. 
The mesh was then refined as shown in Fig. 9a and after a 
further 100 global time steps the density distribution is as 
illustrated in Fig. 9b. The next Level of refinement is 
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displayed in Fig. 10a with the corresponding density contours 
•after a further 100 global steps appearing in Fig. 10b. The 
final refinement is then performed as in Fig. 11a, producing the 
density contours . of Fig. 11b after a further ID global time 
steps. The quality of the final solution is good, showing the 
detached bow shock. 

7. CONCLUSIONS 

Explicit solution techniques for the Euler equations have 
been combined with domain splitting and mesh enrichment 
procedures. Practical examples have been presented which 
confirm the viability of the procedures for the analysis of 
steady problems. For the analysis Df true transient problems, 
the areas requiring refinement will traverse the whole of the 
domain under consideration and the refinement strategy wiLl 
require suitable modification in this case [233 . 
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